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Abstract 

Background: Mathennatical nnodels of dynannical systenns facilitate the connputation of characteristic properties 
that are not accessible experinnentally. In cell biology, two main properties of interest are (1) the time-period a 
protein is accessible to other molecules in a certain state - its half-life - and (2) the time it spends when passing 
through a subsystem - its transit-time. We discuss two approaches to quantify the half-life, present the novel 
method of in silico labeling, and introduce the label half-life and label transit-time. The developed method has been 
motivated by laboratory tracer experiments. To investigate the kinetic properties and behavior of a substance of 
interest, we computationally label this species in order to track it throughout its life cycle. The corresponding 
mathematical model is extended by an additional set of reactions for the labeled species, avoiding any double- 
counting within closed circuits, correcting for the influences of upstream fluxes, and taking into account 
combinatorial multiplicity for complexes or reactions with several reactants or products. A profile likelihood 
approach is used to estimate confidence intervals on the label half-life and transit-time. 

Results: Application to the JAK-STAT signaling pathway in Epo-stimulated BaF3-EpoR cells enabled the calculation 
of the time-dependent label half-life and transit-time of STAT species. The results were robust against parameter 
uncertainties. 

Conclusions: Our approach renders possible the estimation of species and label half-lives and transit-times. It is 
applicable to large non-linear systems and an implementation is provided within the PottersWheel modeling 
framework (http://www.potterswheel.de). 



Background 

Motivation 

An increasing number of biological phenomena are 
described by mathematical models, specifically on the 
basis of biochemical reaction networks [1,2]. The 
dynamic properties of these networks are given by their 
model structure, kinetic parameters, initial values of the 
involved species, and externally specified input func- 
tions. The interpretation of an isolated element of the 
network, e.g. a certain rate constant, has only a limited 
meaning, because its effect can only be understood 
when taking the whole network context into account. 
We therefore seek to introduce two dynamical charac- 
teristics which have a physiological meaning, are 
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intuitive to understand, and capture the system kinetics 
on a higher level of abstraction. The first characteristic, 
the label half-life, applies the half-life concept not to a 
species, but to a virtual label attached to the species. 
The second one, the label transit-time, is the time-per- 
iod it takes for a fraction of labeled entities to pass 
through a subsystem of the network. Both quantities are 
calculated using a novel approach called in silico label- 
ing, which is also introduced in the present work. 

In Silico Labeling and Species vs. Label Half-Life 

In a laboratory tracer experiment, a substance is marked 
to better understand the kinetic properties of the dyna- 
mical system [3]. Different tracer substances have been 
used, e.g. radioactive iodine-125 [4,5] or green fluores- 
cent protein-tagged proteins in combination with fluor- 
escence recovery after photobleaching (FRAP) [6]. A 
good tracer does not hamper the flux of the substance. 
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therefore one can assume that the flux of the tracer 
within a certain reaction is proportional to the flux of 
the original species. This is the key property of the in 
silico labeling approach, where an additional set of reac- 
tions is added to an existing mathematical model 
describing the kinetic behavior of a tracer, called the 
label In contrast to real tracer experiments, the in silico 
method offers the opportunity to define dead-ends, 
avoid double-counting of cycling label, and to restrict 
the label to a sub- network of reactions. This allows ask- 
ing specific questions about the original system, like 
how long it takes for 50% of the molecules of a sub- 
stance to travel along a certain path, while in reality an 
alternative path may exist. In addition, predominant 
paths can be identified in deterministic models as has 
been done previously for stochastic systems [7]. 

Mathematically, the half-life T1/2 of a species is 
defined as the time-period until it reaches half of its 
initial amount assuming no influx. For clarity, we denote 
this time-period as the species half-life (SHL), In non- 
isolated and non-linear processes, this time-period dif- 
fers from the amount of time required for 50% of initi- 
ally existing molecules to be processed. For this, we 
introduce the label half-life (LHL), defined as the half- 
life of the label of a species. Equalities and differences 
between the species and label half-life are displayed in 
Figure 1 and proven in the methods section. 

While for simple systems the species half-life can be 
determined analytically, the symbolic integration of a 



Michaelis-Menten kinetics leads to advanced mathemati- 
cal calculations including the Lambert W function [8]. We 
therefore also provide an automatic and generally applic- 
able numerical method to determine the species half-life. 

Label Transit-Time 

Transit-times are discussed in a variety of fields and 
they are, for example, used to quantify how quickly food 
moves through the gastrointestinal tract [9]. When 
describing the dynamics of Markovian particles, the 
mean transit-time denotes the time spent on average in 
a subsystem [10], while the mean sojourn-time also 
takes into account the probability that the subsystem is 
entered at all [11]. In pharmacokinetics, the so-called 
mean residence time values [12] are estimated based on 
empirical data assuming linear kinetics [13]. Apart from 
linearity, no influx for the species of interest is per- 
mitted. Eventually, the estimation is only applicable to 
observable species. The computation of the mean resi- 
dence time is accomplished by the ratio of the area 
under the first moment curve (AUMC) to the area 
under the curve (AUG) of the concentration-time profile 
of a drug [14]. 

We here introduce the label transit-time (LTT) from a 
source to a target pool in a chemical reaction network 
as the time-period after which 50% of all entities resid- 
ing in the source pool at ^ = 0 have reached the target 
pool at least once. The exact path from source to target 
pool is not important in the unconditioned case. The 



Species vs. Label Half-Life in isolated processes Species vs. Label Half-Life in non-isolated processes 
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Figure 1 Species vs. label half-life. Panel A: The species half-life of a substrate 5 in the reaction S ^ P \s plotted for different reaction types 

(solid lines). Except for processes of order 1, the half-life is time-dependent. Since the substrate is not produced in further reactions, the label 

half-life (dashes) equals the species half-life. Panel B: The substrate 5 participates in a production {A 5) and a processing {S ^ P) reaction. 

Now, the species half-life differs except for a linear processing from the label half-life, because the label flux is proportional to the total flux of 

each reaction and is therefore affected by concentration changes through influx of 5. Both panels: The species half-life has been determined 

analytically and numerically according to the methods section. Matlab scripts to reproduce the plots are available in the additional file 1. 
I ) 
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LTT information could be valuable to estimate the time 
for a drug or an enzyme to reach its site of action. 

Extended Reaction Network 

To determine the label half-life, it is important to distin- 
guish entities residing in the source pool at ^ = 0 from 
other entities entering the source pool at later time- 
points. When calculating transit-times, this discrimina- 
tion has to be applied to all pools and fluxes between 
source and target. To achieve this aim, the species of 
interest is computationally labeled and subsequently 
tracked throughout the dynamical model. The labeling 
is realized by an additional set of reactions describing 
the kinetic behavior of the labeled species, depending on 
the kind of time characteristic LHL or LTT, the source 
species, and potentially a target species. 

In case of label half-life calculations, it is sufficient to 
create labeled reactions for all reactions in which the 
source species is a reactant. In fact, labeled reactions are 
prohibited if the source species is a product; this is to 
avoid double-counting the labeled species. In the case of 
transit-time calculations, for all original reactions in 
which labeled species are involved, a new labeled reaction 
is added. In all labeled reactions with the target species 
being the product, the label is removed and accumulated 
in an artificial pool which is used to determine when 50% 
of the existing label has reached the target. 

The label stays virtually attached to a species through- 
out all modifications of the species, such as phosphory- 
lation or relocalizations, e.g. shuttling into the nucleus. 
While the suggested approach can be implemented in a 
straightforward way for monomeric reaction networks 
with only up to one labeled reactant and product, for 
the general case where the reactions involve multiple 
reactants and products or where labeled species may 
form a polymer, a systematic book-keeping of all possi- 
ble combinations of labeled and unlabeled species is 
required. 

As motivated by laboratory tracer experiments the 
fluxes of the additional system are based on the corre- 
sponding fluxes in the original one, which is explained 
in detail in the methods section. 

Profile Likelihood-based Confidence Intervals 

Recently, we suggested a profile likelihood-based 
approach to determine the confidence intervals on cali- 
brated parameter values in mechanistic mathematical 
models [15]. The same reasoning can be applied in 
order to estimate confidence intervals for the time- 
dependent label half-life and transit-time characteristics. 

Implementation 

All concepts have been implemented within the Potters- 
Wheel modeling and parameter estimation framework 



that is available from http://www.potterswheel.de [16] 
and have recently been applied by the authors to the 
mathematical models of the erythropoietin and epider- 
mal growth factor receptors [17,18]. The application of 
the method within the PottersWheel framework is 
described in additional file 1. 

In the next section, the proposed labeling method is 
illustrated for the JAK-STAT signal transduction path- 
way and afterwards described in detail. After proving 
the equality of species and label half-life for isolated or 
linear processes, a fitted model of the JAK-STAT path- 
way is used to determine the label half-life of unpho- 
sphorylated STAT and its label transit-time when 
cycling through the nucleus of a cell. 

Methods 

Illustration of the method 

Figure 2 illustrates the in silico labeling approach for the 
JAK-STAT signal transduction pathway, where STAT 
molecules cycle between cytoplasm and nucleus. First, 
cytoplasmic STAT molecules {S) are phosphorylated 
{pS) by an active receptor {pR) and form dimers 
{pS_pS), The complexes enter the nucleus {npS_npS) 
where they act as transcription factors, disassociate and 
are dephosphorylated (nS) again. Finally, they return to 
the cytoplasm (S) and can be activated again. In order 
to determine the label half-life of cytoplasmic STAT and 
the label transit-time for a whole cycle, we set source 
and target species to unphosphorylated cytoplasmic 
STAT. At ^ = 0, all molecules of the source pool are 
labeled, symbolized by the small red spheres. The label 
is not removed until the target pool is reached, in this 
case when a STAT molecule leaves the nucleus. Then, 
the label is accumulated in an artificial pool of returned 
label and an unlabeled STAT molecule enters the cyto- 
plasm. Over time, the fraction of labeled to free, unla- 
beled STAT molecules, S^/S^, decreases in the 
cytoplasm. The total flux v]^ of the first reaction, that is 
the phosphorylation of STAT molecules, can be divided 
into the flux v\ of labeled and the flux of free STAT 
molecules. The fraction of the fluxes is set to match the 
fraction of labeled to free STAT molecules, by the fol- 
lowing relationship: 

-| := —Withv[+v{ = vl. (1) 

The label half-life of STAT at time-point t is given by 

LHLs(t) = arg,(^S^{t') = ^)-t. (2) 

The label transit-time from STAT to STAT at time- 
point t can be derived from the time-profile of the 
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Figure 2 In silico labeled JAK-STAT signaling pathway. STAT molecules S (blue) are phosphorylated by an active receptor-kinase complex 
(pR) and form dimers (pS_pS). These dimers enter the nucleus, dissociate, and are subsequently dephosphorylated. Finally, the single STAT 
molecules re-enter the cytoplasm, where they can again be phosphorylated and thus continue the nuclear-cytoplasmic shuttling. The labeling 
approach is visualized by red spheres attached to the STAT molecules. At t = 0, all cytoplasmic STAT molecules are labeled. After the nuclear 
export, the label is removed from the molecule and enters the artificial pool of returned label. Consequently, an increasing fraction of 
cytoplasmic STAT molecules are not labeled which has to be considered in the calculation of fluxes for free and labeled entities. To determine 
the time-dependent label half-life and transit-time values, the labeling procedure is repeated for a series of time-points. 



returned label RL: 

LTTs^sCO = arg, f RL{t') = ^ ) - t. 



(3) 



This procedure is repeated for a series of time-points t 
in order to determine LHL(^) and LTT(^) for all time- 
point of interest. 

Terminology 

In the following we assume that the biological system is 
mathematically described by a set of reactions rp I < J < 
n, corresponding to a set of coupled differential equa- 
tions. The concentration change of each entity Xi, 1 < / 
< m, is the sum over all fluxes of reactions where it 
appears as a product minus the sum over all fluxes of 
reactions where it appears as a reactant, mathematically 
[19] 



= J2 ^ij^i ~ J2 ^ij^i for i < I < m. 



(4) 



Here, Vj describes the flux of reaction j, Uij > 0 the 
stoichiometry of Xi as a product in reaction ; and bij > 0 
the stoichiometry of Xi as a reactant in reaction y. We 



use the same symbol for an entity and its concentration, 
[Xi\ - Xi, The time-profile of each species can then be 
calculated for given initial values = Xi[to) and poten- 
tially driving input functions U/^t), The flux Vj of reac- 
tion J may be a non-linear function of one or more 
species concentrations Xi and externally defined U/^, To 
improve readability, we omit explicitly denoting the 
time-dependency, i.e. Xi{t) is rather written as Xi. 

Analytical and numerical half-life calculation 

The half-life of a species Xi of interest is determined by 
extending the differential equation network (4) by one 
equation for an artificial quantity y depending only on 
the outfluxes of Xi, 



(5) 



with initial value = xj*. The whole system (4, 5) is 
solved either analytically or numerically and the species 
half-life of Xi is given by T1/2 for which 



2 



(6) 
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Note that a half-life characterizes the decay of a quan- 
tity, independent of any production rates. Therefore, all 
influx contributions are neglected in equation (5). In 
general, only linear processes possess a constant half- 
life. Otherwise, the half-life depends on the initial con- 
centration and is therefore time-dependent. In this 
case, the above procedure is repeated for a series of dif- 
ferent initial time-points to. In a numerical integration, 
it is important to limit the maximum integrator step 
size for an accurate approximation of the threshold 
crossing. 

The half-life of a species Xi is only partially related to 
the time it takes for 50% of an experimental tracer to 
leave the source pool. The two values coincide if Xi has 
either no influx or when the outflux from Xi is described 
by a linear process, which will be proved in the next 
two subsections. Therefore, we suggest the in silico 
labeling half-life as a means to determine a time-charac- 
teristic which is motivated by laboratory tracer experi- 
ment with the additional property to avoid tracer- 
double counting in kinetic cycles. 

In silico labeling half-life for isolated processes 

For simplicity, we assume that the species of interest x 
, Xj^ is consumed only in one reaction. In in 
silico labeling, the flux of the corresponding label z 
depends on the outflux of x by 



(7) 



The in silico labeling half-life of x is defined as the 
time when z drops to Zo/2, We will show that this time 
equals the species half-life of x if its influx is zero. 
This property is independent from the amount of initi- 
ally labeled entities, i.e. it holds for any Zq/xq g R"^: 

Proof: 

Let X be determined by the processing with an 
unknown, potentially non-linear outflux Vout and no 
influx Viyi = 0, i.e. v = v^^^, 

x = —V and Xq = x[0). 

Then, the kinetics of the label species z(t) is given by 
z 

z = — v[x) and zq = z[0). (8) 

X 

It can be shown that the factor § is constant: 



dt 



zvx + zv 



Since this relation holds also true for ^ = 0, the pro- 
portionality constant is given by / = ^- Then, equation 
(8) reads 

z(t) = —fv and zq = fxQ. 

Both processes x{t) and z{t) share the same half-life 
ri/2, since 

^(Ti/2) = y ^ 4^1/2) =MTi/2) =/y = y , q-e.d. 

This relation does not hold for processes with ^ 0, 
because the fraction z/x becomes time-dependent as the 
labeling gets diluted, except for linear outfluxes as 
shown in the next section. 

In silico labeling for linear processes 

In this section, we prove that the label half-life coincides 
with the half-life of a species x which is produced by an 
unknown, potentially non-linear influx V/„ and is con- 
sumed by a linear process. 
Proof: 

Let X be given by an unknown, potentially non-linear 
influx Vjyi and a linear outflux, kx, 



kx 



Then, the analytical half-life of x can be determined 
via 



y 

Tl/2(X) 



—ky 
ln(2) 



For the labeled system z it holds that 



-kx 



X 

-kz 



=^z[t) 



kt 



ZqC 

ln(2) 

k 



q.e.d 



Creating the Extended Reaction Network 

Some entities Xi belong to the group of tracked, i.e. 
potentially labeled entities. Let us assume that they are 
given by jvi, . . . , Xa and untracked ones by Xa+i> . . . , 
x^. Further, it can be assumed without loss of generality 



that (1) Xi, 



Xy < a are not complexes consisting of 



two or more tracked single entities, and (2) that the 

tracked single entities within each complex 

belong to the set . . . , Xy, In the JAK-STAT example. 
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S, pR_S, and pS belong to Xi, , , , y Xa and pS _jfS to Xa 
+1,. . . , Xfn^s it contains two labeled single entities pS. 

Creating additional entities x^'' 

A new set of labeled or free entities x^^ is created based 
on the original x, by applying the following rules: 

♦ Start with an empty set, x^^ = {} 

♦ Single entities: For each Xi g {xi, . . . , Xy}, x^^ is 
enlarged by a labeled x[ and a free xf version of the 
original entity 

♦ Complex entities: Each complex Xi g {Xy+i, . . . , 
Xa} is decomposed into n\xi, . . . ,nyXy. Due to the 
combinatorial multiplicity, 

2^H";- (9) 



possible combinations using labeled x^^ and free x- 
entities are created, taking into account the order of 
the elements in the original complex Xi, and are 
added to x^^. The complex pS pS for instance leads 
to the four new complexes pS^_pS^, pS^_pS^, 
pS^_pS^, and pS^_pS^, 

Creating additional reactions r^'' 

In order to create a new set of reactions r^^ , the com- 
binatorial multiplicity has to be applied not only to 
complexes but also to the ordered lists of reactants 
and products. Suppose an ordered list / of entities 
from the set {Xi}i < with possible repetition, as for 
example the reactants of the reaction A + A + pA_pA 
A_A_pA_pA corresponding to / = {A, A, pA_pA), 
Summing up all single reactants and elements of the 
complexes leads to p single entities, in this case = 4. 
Taking into account all combinations of labeled and 
free entities, 2^ different lists can be derived, in the 
example 

h = (A^, A^, pA^jpA^), h = (A^ A^ pA^jpA^), and lie = (A^ A\pA^.pA^). 

Without loss of generality, only the first S reactions of 
the original system are assumed to affect a tracked 
entity. In these reactions, at least one reactant or pro- 
duct is a tracked entity. Then, a new set of reactions r^^ 
can be established. Starting with the empty set r^^ = {}, 
for each reaction r/ g {ri, . . . , r^} with one or more 
reactants of tracked entities, 

1. all reactants and products not belonging to the 
group of tracked entities are removed. 



2. the combinatorial multiplicity approach is applied 
to the ordered list / of the remaining reactants lead- 
ing to Ii, ...J2P> 

3. 2^ reactions are added to r^^ with reactants 
Ii, ...,I2P and the corresponding products. 

4. the fluxes v{^,...,ij^ of the new reactions are 
given by 

^ n^eif^ ^ ^ ^^^^ 



Note that again the same symbol has been used for 
the entity name and its concentration. The sum over 
all weighting factors is 1. 

Reactions g {ri, . . . , r^}without reactants produce 
only free entities, which simplifies the conversion of 
before adding to r^^: All untracked entities are removed, 
all Xi are replaced by xf , and the flux is again given by 
equation (10). 

When calculating the label half-life, products that 
coincide with the initially labeled entity are replaced by 
the corresponding free entity. This corresponds to 
removing the label and is necessary to avoid double- 
counting and to exclude upstream fluxes. 

In order to calculate the label transit-time, entities 
entering the target pool must be released from their 
labeling, again, to avoid double-counting. Therefore, all 
labeled target entities are replaced in the reaction net- 
work r^^ by their free counterparts. At the same time, a 
new product is added to those reactions where the tar- 
get entity is a product to accumulate the returned label, 
RL. 

Calculating the Label Half-Life and Transit-Time 

Since the label half-life and transit-time characteristics 
are time-dependent, the label is not only injected at 
time-point 0, but the procedure is repeated for a series 
of time-points t (let Xi be the source species): 

1. Set all initial values for labeled entities and RL, if 
available, to 0. Set the initial value of free entities to 
the value of their counterpart in the original 
network. 

2. Numerically integrate the ordinary differential 
equations corresponding to the extended reaction 
network {r, r^^} from 0 to t, 

3. Apply a complete labeling of the source species: 
Set x[(t)=Xi(t) and xf (t) = 0. This step corre- 
sponds to the label injection. 
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4. Continue the numerical integration. 

Threshold crossing at t" of the time-profiles x^[f > t) 
and RL(^' > t) with x^[t)/2 defines the label half-life and 
label transit-time as t"'t, respectively. The threshold 
crossing is determined by linear interpolation of the dis- 
crete samples given by the numerical integration. 

Profile Likelihood-based Confidence Intervals 

We recently suggested a profile likelihood-based 
approach to determine simultaneous and separate confi- 
dence intervals for calibrated unknown model para- 
meters [15]. In order to determine confidence intervals 
for the calculated label half-life and transit-times, the 
above procedure is not only repeated for a series of 
time-points, but also for a series of parameter settings. 
Each setting corresponds to one extreme point on the 
multi-dimensional manifold of acceptable parameter 
values, where one parameter has reached a lower or 
upper confidence threshold. By plotting all LHL or LTT 
profiles into one axis and creating an envelope between 
the largest and lowest values, a confidence interval for 
LHL and LTT is given. 



In order to calculate the half-life for Michaelis-Menten 
kinetics, x = —Vmax^/iKm + x), the following integral 
form is used which has been derived in [20], for x{t) 
with known Xq dX t = to: 

Xo-X + Km\n (^^^ = Vmax[t to) (17) 

Panel A of Figure 1 displays the analytic results and 
their numerical approximation. 

Results 

In this section, the in silico labeling approach is applied 
to the JAK-STAT signaling pathway. The following 
mass action-based mechanistic model of the pathway 
has been calibrated to immunoblot measurements for 
Epo-stimulated BaF3-EpoR cells (model motivated by 
and data taken from [21]): 

d{^)ldt = -feiSpR + fesnS (18) 
d[pS)/dt = feiSpR-2fe2pSpS (19) 
d{pSjpS)/dt = kj pSpS — fespS-pS (20) 



Analytic half-lives for simple, isolated processes 

The half-life Ti/2{t) of simple and isolated biochemical 
reactions can be calculated analytically. Except for first- 
order processes, it usually depends on the concentration 
^0 = ^{to) at the time-point of interest to and is therefore 
time-dependent: 

Process of order 0: T1/2 = ^ (11) 



Process of order 1: T1/2 = — (12) 



Process of order 2 : T1/2 = (13) 

Xo -k 

2"-i - 1 

Process of order n > 1 : T1/2 = — r (14) 

[n - l)k -Xq 

Michaelis - Menten : T1/2 = (15) 

^max 

The half-life calculation for a process of order n >1 
with X = —kx^ is based on the integral form 

= 3Pr + (" - 1)'^^- (16) 

A Xq 



d(npS_npS)/dt = fespS-pS — fe4npS_npS (21) 

d[nS)/dt = 2k4 npS_npS - ks nS (22) 

A smoothing spline approximation of the phosphory- 
lated receptor served as the input function pR{t) trigger- 
ing the phosphorylation of STAT {S pS), After 
dimerization {pS + pS ^ pS.pS), the complexes enter the 
nucleus {pS_pS npS_npS), Then they dissociate and 
are dephosphorylated {npS_npS nS + nS), Finally, sin- 
gle STAT molecules leave the nucleus again {nS S), 
Model parameters were estimated using a Levenberg- 
Marquardt approach and the PottersWheel modeling 
software. The pools of total and phosphorylated cytoplas- 
mic STAT have been used as observation functions. The 
kinetic parameters were estimated as /q = 1.37, /c2 = 0.22, 
/C3 = 0.63, /C4 = 0.59, and ks = 0.59. The initial value of S 
was calibrated to 0.96 and the scaling factors for the 
observables to 1.45 for pS.obs and 0.98 for S.obs. 

Labeled system 

In order to determine the label half-life and transit-time 
of STAT, S is both, the initially labeled entity and the 
target pool. The flux of the label is illustrated in Figure 
2. The time-courses of the original (solid blue) and 
labeled system (dashed red) are compared in Figure 3. 
In the beginning, both systems behave in the same man- 
ner. Then, the first wave of STAT molecules return 
from their cycle through the nucleus. Since they loose 
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nS Initial and returning label 




Time(min) 



Figure 3 Time-courses of the original and labeled system. At t = 0, all STAT molecules are labeled and the original system (solid blue) 
coincides with the labeled one (dashed red). After the initial signal wave, STAT molecules return from the nucleus to the cytoplasm. Since the 
label transit-time for a complete cycle of cytoplasmic STAT is investigated, the returning molecules release their label, which is counted in an 
artificial entity (bottom right, green). Unlabeled molecules are depicted in yellow. Solid rose lines depict half labeled species, i.e. dimers of a 
labeled and an unlabeled molecule. They are shown only once, since the trajectory for example of pS^.pS^ is identical to the one of pS^.pS^. 



their label, the amount of labeled cytoplasmic STAT 
does not recover in contrast to the amount of STAT. 
After 13 minutes, 50% of the initially labeled STAT 
molecules passed the nucleus at least once, as shown by 
the artificial pool of the returned label. The bimodal 
behavior of pSTAT exemplifies the first original signal 
wave and the secondary cycling effects. The in silico 
labeling approach allowed for discrimination between 
these two dynamics. In order to determine the transit- 
time for t >0, the label is injected at a series of time- 
points, which is visualized in Figure 4. 

Label half-life and transit-time 

Figure 4 depicts the label half-life of STAT as calculated 
from the time-course of S^, It reaches a minimum of 0.6 



±0.1 minutes after ten minutes compared to a half-life 
of approximately 3 to 4 minutes at the starting point of 
the time-course analysis. For later time-points, the sti- 
mulus decreases (not shown) and STAT is no longer 
phosphorylated, resulting in an increased label half-life 
of STAT molecules. The minimum label transit-time for 
a complete cycle of STAT molecules was estimated as 
12 ± 2 minutes. 

Profile likelihood-based confidence intervals 

In order to investigate how uncertainties in calibrated 
model parameter values propagate to the estimated time- 
characteristics, we applied the profile likelihood approach 
on an identifiable version of the model. The kinetic para- 
meters involved in phosphorylation (/ci), dimerization {k'2), 
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Trajectories for labeled STAT 
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Figure 4 Time-dependent label half-life and transit-time of cytoplasmic unphosphorylated STAT. A: In order to determine the label half- 
life and transit-time, a family of labeled STAT trajectories is calculated. For trajectory /', the amount of S(t/) label is injected into the system as 5^ 
at ti, with 0 < ti < 30. Therefore, the upper limit of subplot A matches the time-course of S for t < 30 in Fig 3. B: Shown are the corresponding 
trajectories for returned label. The label half-life for trajectory / is the time period until drops below S^[ti)/2. C: The shortest label half- 

life is given as 0.6 ± 0.1 minutes. D: The label transit-time for a complete cycle of STAT molecules through the nucleus is the time-period until 
the returning label exceeds half of its initial value. Here, the minimum was estimated as 12 ± 2 minutes. Previously, it has been shown that the 
sojourn time of a single STATS molecule in the nucleus is about 6 minutes [21]. Our results indicate that on (median) average, a STATS molecule 
spends equal times in each compartment and requires about 12 min for one cycle from the cytoplasm to the nucleus and back. The 9S% 
confidence intervals have been determined using the profile likelihood approach (RLE) and are displayed as grey envelope curves. 



nuclear import (k^) and export {ks) were systematically 
varied consecutively within four orders of magnitude 
between O.Olkfit and lOOkfn, with kfn being the parameter 
value for the best fit. For each variation, the other free 
parameters were calibrated resulting in a profile likelihood 
estimation (see Fig. S2 in additional file 1). All parameter 
settings corresponding to a crossing of the profile likeli- 
hood with the X^-threshold of the separate 95% confi- 
dence interval are used to recalculate the label half-life 
and transit-time. Figure 4C and 4D display the LHL and 
LTT 95% confidence interval by envelope curves. In case 



of the label half-life of cytoplasmic STAT, the confidence 
interval is very narrow allowing the LHL estimation within 
±0.1 minute for a range of label injection times between t 
= 0 and t = 20 minutes. The label transit-time has a wider 
confidence interval reflecting the larger number of reac- 
tions involved in a complete cycle of shuttling STAT. 

Discussion and Conclusions 

In this paper, the half-life of a species has been com- 
pared conceptually, analytically, and numerically to the 
half-life of a label in a hypothetical tracer experiment. 
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Two time-characteristics, the label half-life and label 
transit time have been introduced, which capture the 
kinetics of a dynamical system on a higher level than e. 
g. single rate constants. Calculation of the time-charac- 
teristics and their profile likelihood-based confidence 
intervals for an identifiable pathway model showed that 
the approach is robust against parameter uncertainties. 
The quantities are calculated based on the novel in silico 
labeling method, which relies on an extended reaction 
network taking into account constraints concerning 
double-counting, upstream fluxes and combinatorial 
multiplicity. Our model-based in silico approach allows 
for insights into reaction networks that cannot be deter- 
mined experimentally. 

The proposed method provides important information 
for a wide spectrum of biological applications ranging 
from cell biology and pharmacokinetics to population 
dynamics. We applied it to a non-linear model of the 
cellular JAK-STAT signaling pathway, which allowed for 
calculating the time-dependent label half-life and tran- 
sit-time of cytoplasmic STAT. 

In summary our approach enables to calculate the 
amount of time a molecule spends in a certain state or 
compartment and therefore provides novel insights into 
the temporal scale of networks. This knowledge will 
have profound impact on drug design, as it offers the 
possibility to predict the life-time of a specific molecule 
and provides a basis to improve drug targeting. 

Additional material 



Additional file 1: Application within PottersWheel. This additional file 
contains MATLAB scripts to run various tasks related to the in silico 
labeling approach, http://www.biomedcentral.com/imedia/ 
4654854926777309/suppl .pdf.. 
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